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, Based on the statistical dynamical mean field theory, we investigate, in a generic model for a 

JlJ . strongly coupled disordered electron-phonon system, the competition between polaron formation and 

' Anderson localization. The statistical dynamical mean field approximation maps the lattice problem 

to an ensemble of self-consistently embedded impurity problems. It is a probabilistic approach, 
f — . focusing on the distribution instead of the average values for observables of interest. We solve the self- 

CNj ' consistent equations of the theory with a Monte-Carlo sampling technique, representing distributions 

for random variables by random samples, and discuss various ways to determine mobility edges from 
the random sample for the local Green function. Specifically, we give, as a function of the 'polaron 
D ' parameters', such as adiabaticity and electron-phonon coupling constants, a detailed discussion of 

, ' the localization properties of a single polaron, using a bare electron as a reference system. 
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PACS numbers: 



I. INTRODUCTION 



^ . . . . . 

^ ' There is a large variety of materials in which due to strong electron-phonon coupling electrons and phonons lose 

O their identity and form new composite entities: polarons. Examples of current interest are, among others, the high- 
I ^ I temperature superconducting perovskites (Bar- Yam et al. 1992, Salje et al. 1995), the non-metallic nickelates and 
] bismuthates (Cheong et al. 1994, Oyanagi et al. 2001), the colossal magneto-resistance manganites (Jin et al. 1994, 
. Tokura et al. 1999), and semiconducting molecular crystals (Parris et al. 2001). In all these materials, electronic 
^ ' properties, for instance the electronic spectral function or the optical conductivity, show substantial polaronic effects. 
, The most important signature in this respect is thermally activated transport (along at least one of the crystallographic 
1^*^ axes j . 

I , The polaron concept has been introduced by Landau (1933) and ever since played a central role in the analysis of 
strongly coupled electron-phonon systems. In a deformable lattice, the coupling between an electron and the lattice 
leads to a lattice deformation whose potential tends to bind the electron. This process is called self-trapping, because 
■ the potential depends on the state of the electron (Firsov 1975, Rashba 1982, Wellein and Fehske 1998). Translational 
symmetry is not broken in this process, transport at low enough temperatures is still coherent and band-like, with 
a much larger mass, however, because of the lattice distortion the electrons have to carry along. The self-trapping 
process, i.e. the formation of polarons, is therefore only the crossover from a weakly dressed electronic quasi-particle to 
a heavily dressed polaronic quasi-particle. Even the small polaron, although much less mobile than the bare electron, 
is still itinerant and not localized. 



I 

Q _ Details of the actual materials notwithstanding, the Holstein model (Holstein 1959), which describes electrons locally 
O ■ coupled to dispersionless phonons, captures the essence of the self-trapping process. Most of the theoretical work is 
^ ' therefore directly based on the Holstein model (DeRaedt and Lagendijk 1983, Ranninger and Thibblin 1992, Marsiglio 
• ^ 1993, Alexandrov et al. 1994, Stephan 1996, Wellein et al. 1996, Capone et al. 97, Bonca et al. 1999). Extensions 
to long-range electron-phonon couplings on a discrete lattice were discussed by Alexandrov and Kornilovitch (1999) 
and by Fehske et al. (2000). Despite the extensive numerical simulations, the properties of the Holstein model are 
not yet fully understood, especially at finite densities, where the situation changes because of the interaction between 
polarons. But even in the extreme dilute limit, where polarons in first approximation do not interact, many questions 
remain. For instance, little is quantitatively known about polaron formation in a disordered environment. 

On the other hand, polaronic materials are complex materials, where chemical as well as crystallographic imper- 
fections can be quite substantial. Accordingly, the self-trapping process, producing a heavy quasi-particle, occurs in 
an environment where, strictly speaking, translational symmetry is broken, and where defects most probably act as 
nucleation sites for the formation of polarons. Therefore, a complete theory of polaron formation has to take the 
disordered environment into account. In particular, the notoriously difficult question of whether polaron states are 
band states or localized defect states can be only meaningfully addressed in a context which takes the most important 
effect due to disorder, the possibility of Anderson localization, explicitly into account. 

In his seminal work Anderson (1958) has shown that disorder, if sufficiently strong, dramatically changes the 
properties of the electronic states in a solid. Whereas for small disorder the electronic states are extended Bloch waves, 
for large disorder electronic states are localized defect states. Anderson localization, i.e., the transition from extended 
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to localized states, is well understood for noninteracting electrons [see, for instance, the review articles by Thouless 
(1974), Brczini and Zcrki (1992) and Kramer and MacKinnon (1993)]. It depends strongly on dimensionality. In one 
dimension [and arguably in two dimensions (Abrahams et al. 1979)] localized and extended states cannot coexist, 
because infinitesimally small disorder suffices to localize all states at once, whereas in three dimensions mobility 
edges separate localized from extended states. Since the mobility edge is of utmost importance for the transport 
properties (Mott 1968a, Mott 1968b, Mott 1976, Mott 1981), an enormous amount of experimental and theoretical 
effort has been directed towards a precise determination of its position. 

That Anderson localization might affect polaron formation, and vice versa, has been emphasized several times, 
starting with Anderson (1972) himself, who pointed out that the mobility edge might be surrounded by self-trapped 
states. Not much quantitative work exists however to determine, for instance, the position of the mobility edge and 
its consequences for the conductivity (Girvin and Jonson 1980, Mueller and Thomas 1983), or the character of the 
states near the mobility edge (Cohen et al. 1983) taking electron-phonon coupling explicitly into account. Despite 
some occasionally fruitful adaptations of techniques from nonlinear dynamics (Kopidakis et al. 1996), where polaron 
formation is considered as an intrinsic interaction-driven nonlincarity, the precise mechanisms of the interplay of 
Anderson localization and polaron formation are essentially unknown. 

As a first step to address this problem with techniques developed in the field of strongly correlated electron systems, 
Bronold et al. (2001) applied the dynamical mean field theory (DMFT) (Georges et al. 1996) to study a single electron 
in the Holstein model with binary disorder, focusing especially on polaron states at the high energy edge; of the lowest 
polaronic subband. These states are extremely sensitive to disorder, with a return probability which, due to the 
phonon admixture, is several orders of magnitude larger than for the states at the bottom of the subband. Within the 
DMFT it was however impossible to determine the critical disorder needed to localize polaron states. Accordingly, 
mobility edges for polaron states could not be determined. 

To overcome the limitations of the DMFT, we adopted therefore a recently proposed generalization of the DMFT, the 
statistical dynamical mean-field theory (statDMFT) (Dobrosavljevic and Kotliar 1997, Dobrosavljevic and Kotliar 
1998), to the Holstein model with uniformly distributed on-site energies (Bronold and Fehske 2002, Bronold and 
Fehske 2003). The statDMFT is not only capable of accounting the polaron formation process but also the spatial 
fluctuations giving rise to Anderson localization. Within the statDMFT we could clearly distinguish between itinerant 
and localized polaron states. The purpose of this paper is to give an extended account of our investigations, with 
a complete description of the technical apparatus, including the details of the numerical implementation of the 
statDMFT, and a full discussion of the localization properties of a single polaron. 

As a generic model for a disordered polaronic material, we use the Anderson-Holstein model (AHM) 

H = Y^ eiJiia - Jijclcja + ^Y1 " \/%^Yibi + bl)ni„, (1) 

where c|^(6|) are electron (phonon) creation operators, rijo- = cj^Cjo- is the electron density on site i, the electron 
transfer integral Jij = J for {i, j) next neighbour sites and zero otherwise, Cl is the bare phonon energy {h = 1), 
and Ep is the polaron shift. The on-site energies {ei} arc assumed to be independent, identically distributed random 
variables with a uniform distribution p{ei) = {l/^)9{"f /2 — \ei\). 

The 'polaron properties' of the AHM are governed by two parameter ratios: the adiabaticity a = il/J, indicating 
whether the polaron is light a ^ 1 or heavy a » 1 (Rashba 1982, Wellein and Fehske 1998), and a dimensionless 
electron-phonon coupling constant, A = Ep/2J or = Ep/Q.. Polaron formation sets in if both A > 1 and > 1. 
The internal structure of the polaronic quasi-particle depends on a, in particular the momentum dependent phonon 
admixture. It is quite different in the adiabatic (a <C 1), non-adiabatic (a ~ 1), and anti-adiabatic (a ^ 1) regimes. 
Another important effect occurs in the intermediate coupling regime, where electron-phonon coupling initiates long- 
range timncling processes (Fehske et al. 1997a, Fehske et al. 1997b). Evidently, how disorder affects polaron states 
strongly depends on the polaron parameters. As a result, the 'localization properties' are expected to be quite different 
in the various polaronic regimes. 

The strength of disorder is specified by the width 7 of the distribution for the on-site energy e^. The structure of 
the AHM suggests to distinguish two regimes: The weakly disordered Holstein regime, where 7 is small on the scale 
of the bare bandwidth, and the strongly disordered Anderson regime, with 7 large on the scale of the width of the 
polaronic subbands. 

The organization of the paper is as follows. Section II gives a general description of the statDMFT. As a guidance, we 
first put in Subsection II. A our work and the statDMFT into perspective. For electrons interacting with phonons, the 
basic equations are derived in Subsection II. B. The statDMFT reduces for lattices with large coordination number to 
the DMFT. This is shown in the appendix. For practical calculations it is convenient to use a Bethe lattice. Subsection 

II. C specifics therefore the basic equations for the Bethe lattice. Subsection II. D finally introduces the localization 
criterion we used to identify localized states. Numerical results in the single particle sector are presented in Section 
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III. After specializing the basic equations in Subsection III. A to the case of a single particle at zero temperature, we 

first give in Subsection III.B results with electron-phonon coupling turned off. The pure Anderson case serves two 
purposes: First, it demonstrates, for a well understood case, that the method and the localization criterion work; 
second, it is used as a reference point for the analysis of the polaron case presented in Subsection III.C. Finally, we 
conclude in Section IV with a short outlook. 



II. STATISTICAL DYNAMICAL MEAN FIELD THEORY 

A. General considerations 

Before we go into technical details, we put our work and the statDMFT into perspective. Self-trapping in a 
disordered environment is one example of physical problems where interaction and disorder effects interfere with each 
other. The most prominent example is electron transport in the impurity band of doped semiconductors. Especially 
near the metal-insulator transition, Coulomb interaction and disorder effects cannot be separated, and it is suspected 

that it is precisely the interplay of both, which leads to new physical phenomena, for instance, to the formation of 
local moments (Finkel'stein 1983, Castellani et al. 1984) or to a new emerging magnetic phase (Kirkpatrick and Belitz 
1990). 

There is no method which could be mechanically used to investigate disordered interacting systems beyond pertur- 
bation theory. In particular, the strongly disordered regime, where Anderson localization takes place, is beyond the 
applicability of most theoretical techniques. Nevertheless, there have been various attempts to generalize methods, 
which have been successfully employed to investigate Anderson localization of noninteracting electrons, in such a way 
that they can be also utilized in situations where interactions are important. Most efforts have been directed towards 
electrons interacting via short- and/or long-range Coulomb potentials, with a few exceptions dealing with electrons 
coupled to phonons (Anderson 1972, Girvin and Jonson 1980, Cohen et al. 1983, Mueller and Thomas 1983, Kopidakis 
et al. 1996, Bronold et al. 2001, Bronold and Fehske 2002, Bronold and Fehske 2003). 

After Anderson's original work (1958) and Mott's suggestion (1968a, 1968b, 1976, 1981) of its connection to trans- 
port properties in amorphous semiconductors, three major techniques have been cultivated to study Anderson localiza- 
tion quantitatively in further detail: (i) direct numerical simulations [see, e.g., the review articles by Thouless (1974), 
Brezini and Zerki (1992), Kramer and MacKinnon (1993)], (ii) diagrammatic techniques [VoUhardt and Wolfle (1980), 
for recent studies consult also the review articles by Kramer and MacKinnon (1993) and VoUhardt and Wolfle (1992)], 
and (iii) field-theoretical approaches [Wegner (1976), for recent studies consult Kramer and MacKinnon (1993)]. Sub- 
sequently, all three techniques have been applied to interacting systems as well. Direct numerical simulations are 
conceptually the simplest approach. However, they are restricted to small system sizes and usually to a truncated 
configuration type treatment of the Coulomb interaction. To this c;atc!gory belong, for instance, all attempts to study 
the localization of two interacting electrons in the background of a frozen Fermi sea (von Oppen et al. 1996). In the 
case of electron phonon interaction, the size of the phonon Hilbert space is the limiting factor. Diagrammatic methods 
work usually only for weak interactions (Lee and Ramakrishnan 1985). In principle, strong interactions can be also 
treated with diagrammatic techniques, if they are based on renormalized expansions, but the proliferation of diagrams 
makes the practical implementation very often rather complicated if not impossible. Field-theoretical techniques, in 
contrast, for instance based on the construction of an effective field theory for the low-energy, long- wavelength exci- 
tations (Belitz and Kirkpatrick 1994, Kamenev and Andreev 1999), albeit very promising, operate in abstract spaces, 
with order parameter functions whose physical content is sometimes hard to grasp. 

Anderson (1958), on the other hand, used a completely different approach. Instead of focusing on the calculation 
of averaged correlation functions, he emphasized that, in general, average values are not representative, especially 
in the strongly disordered regime. All variables of the theory, in particular, observables, should be characterized 
by distributions. Specifically, he examined, in probabilistic terms, the convergence properties of the renormalized 
perturbation series for the local hybridization function appearing in a locator expansion of the local Green function. 
Towards that end, he gave asymptotic estimates for the higher order terms of the series. The mathematical arguments 
are quite involved and have been further clarified later, most notably, by Thouless (1970, 1974) and Economou 
and Cohen (1972). A simplified probabilistic approach, based on the self-consistent solution of the second order 
renormalized perturbation series, started with Abou-Chacra, Anderson and Thouless (1973). It proved to be a very 
powerful tool, applicable not only to substitutionally (Kumar et al. 1975, Brezini 1982, Brezini and Olivier 1983, 
Miller and Derrida 1993) but also to topologically disordered materials (Heinrichs 1977, Fleishman and Stein 1979, 
Elyutin 1979, Elyutin 1981, Logan and Wolynes 1985, Logan and Wolynes 1986, Logan and Wolynes 1987). 

The flexibility of the simplifled probabilistic analysis of the renormalized perturbation series has been recently 
utilized by Dobrosavljevic and Kotliar (1997, 1998) to analyse various physical properties of disordered electrons with 
strong local electron-electron correlations (Miranda and Dobrosavljevic 2001a, Miranda and Dobrosavljevic 2001b, 
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Aguiar et al. 2003). They successfully combined the dynamical mean field approximation for the description of (local) 

correlations with Abou-Chacra, Anderson, and Thoulcss's probabilistic treatment of the spatial fluctuations due to 
randomness. The statDMFT maps the original lattice problem onto an ensemble of impurity problems, which is then 
analysed in terms of probability. In contrast to the field-theoretical approaches, the method is very intuitive and 
conceptually simple. Moreover, it captures local interaction processes nonperturbatively. It is therefore capable of 
addressing issues beyond standard diagrammatic perturbation theories. Numerically it is somewhat involved, but not 
restricted to small system sizes as direct muncirical simulations. The main drawback of the method is its restriction to 
local, short-range interactions. In principle, long-range interaction processes could be incorporated, but only at the 
level of static molecular fields (for instance, in the case of Coulomb interaction a Hartree term would appear), which is 
probably not sufficient for a complete investigation of the interplay between long-range interactions and disorder. For 
the AHM, however, where the (electron-phonon) coupling is local, the statDMFT's problem to deal with long-range 
interactions does not occur and the statDMFT is expected to capture the essential physics of the interplay between 
Anderson localization and self-trapping. 

B. Derivation of the basic equations 

The derivation of the statDMFT for a generic model of interacting electrons has been given by Dobrosavljevic and 
Kotliar (1997, 1998). It is based on the cavity method and applies to arbitrary lattices, temperatures, and densities. 
To make our presentation self-contained and to fix our notation, we outline in this subsection the cavity method, 
carefully paying attention to the probabilistic interpretation of the statDMFT and to specifics due to electron-phonon 
coupling. 

The starting point is the partition function 

Z= j l[VclVci,VblVbie-''^'''''''''''\ (2) 

where the action, for a fixed reahzation of the on-site energies, is given by 

S[c, c\ b, 6t] = I^tJ^ [clir)[dr - ^AcjAr) + bl{r)drbj{T) + H[c, c^, b, 6t](T)| . (3) 

All sites but one are now integrated out, yielding an effective single-site action. For that purpose, the action is split 
into three parts: S = Si + 5^'^ -|- A5(i), where Si is the action of the isolated site i, 5^'^ is the action for the system 
with site i excluded, and 

ASH) = ~ f^drY, [cl{r)riUr) + r?l(r)c;.(T)|, (4) 

with rjia = JiiCia, is the part of the action which connects site i to its neighbours. 

Identifying AS{i) as a source term to generate cavity Green functions, i.e. Green functions for the lattice with site 
i removed, the effective single-site action S{i), defined by 

= i / n P4Pc,-,P6]P6,e-[^-+^"'+^^(^)], (5) 

can be written as 

S{i) = S^ + S, [cl , , b\ , - Vr(^) [r?t , ,7] , (6) 

where we explicitly indicated the fields on which the various terms depend. The first term 5*0 is a constant. While 
the second term, the action for the isolated site i, contains only the fields on site i, the third term, 

W^%\n\ = l+E r^^l- f f dT[... f dT'^A^^-\Tr,...,T'^), (7) 

„ Jq Jo Jo Jo 

with 

AP"l(ri,...,r;) = (-)" Y ■ E '?L.(n)--4.aJ^")G[l...-.;(ri,...,r;)r?,,.;(r0...r7z;.;(r;), (8) 
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involves fields from all sites. 

The effective single-site action S{i) contains the full information of the original lattice action. Integrating out all 
sites but one introduces a 2n— point cavity Green function g|*|^^ ^, (ti, r^). This function is extremely involved 
and approximations are required to obtain feasible equations. A crucial simplification can be made by truncating Eq. 
O after the n = 1 term. This leads to a quadratic effective single-site action, similar to the DMFT (Georges et al. 
1996), where it was further shown that the truncation becomes exact for lattices with infinite coordination number. As 
a consequence, although the cavity construction is performed for lattices with a finite coordination number, interaction 
processes, here due to the electron-phonon coupling, are treated as if the lattice had infinite coordination number. As 
a result of the truncation, the statDMFT becomes a mean field theory. 

Keeping only the n—1 term, Eq. Q depends only on the electronic fields on site i and reduces to 



dT 



dr' V c, 



i(r)i?,<Tm(T- 



(9) 



with the hybridization function 



Hi, 



2 Y^^(i) 

/ J la-ma 
Im 



{r-r'). 



(10) 



The information of the lattice is now contained in the two-point cavity Green function G\]^^^{t — t'), describing all 
paths from site / to site m without passing through site i (note that due to the definition of Jim, both sites / and m 
are neighbouring sites to i). 

Absorbing irrelevant constants, which do not affect the electron dynamics, into Zo, the effective single-site action 
reads 

S{i) = ^^dTdT'^cL(r)|[9, + e,-^]<5(T-r') + ff,;a..(r-r')|c..(r') 

+ + (11) 

with 



(12) 



and 



S^nt{i) = - dT^/E;j}[h{T) + foj (r)]n,, (t) . 



(13) 

According to Eq. Hll|l . the interaction between electrons and phonons takes place only on site i, which, due to the 
hybridization function Hidair), is however embedded into the full lattice. 

Although, for a given Hi„ia{T), any local Green function could be obtained from S(i), the set of equations is not 
closed. This can be most clearly seen from the identity (Dobrosavljevic and Kotliar 1997, Dobrosavljevic and Kotliar 
1998) 

dT'Gl,,,{T - t')G,,„,(t') = j^^ dT'\Gl,,ma{T - t') - G«„,(t - t') | G,,,, (t') , (14) 

which connects the cavity Green function g|^''jjj^(t) with the lattice Green function Gia-ma-iT). Therefore, the hy- 
bridization function Hiaia{T) defined in Eq. (|10|l is a functional of the full, nonlocal lattice Green function, which of 
course cannot be obtained from the effective single-site action S{i). 

To close the set of equations, Dobrosavljevic and Kotliar (1997, 1998) suggested to approximate in Eqs. H1U|I and 
(|14|l . for a given realization of disorder, the full lattice Green function Gia-m<j{T) by the bare lattice Green function 
G'lamrr J with on-sitc energies shifted by the local self-energy. In Matsubara space the substitution reads 



G,^ja{iu}n) G^ „(w„) 



(15) 



FIG. 1: K = 2 Bethe lattice. 



with 'Eia-{iLu) defined by 



lUJn - ei+ fi- Hi„ia[lU}n) - S]icr(jW„) 

which is the Fourier transformed Dyson equation for the local Green function Giaia{T) = —■{TTCia{T)c\^{{)))g^i-^ as 
obtained from the effective single-site action S{i). Note, the substitution Eq. (|15|l is not exact, because in general the 
interaction self-energy is nonlocal. Besides the truncation of the effective single-site action, this is the second major 
approximation invoked by the statDMFT. 

Now the set of Eqs. pU|) ~ (|16|l is closed and constitutes the basic statDMFT equations, applicable to arbitrary 
temperatures, densities, and lattices. A self-consistent solution yields the local Green function Giaia{i^n)- 

It is crucial to realize that, due to the randomness of the on-site energies e^, the local Green function Giai<j{T), 
the local self-energy I]io-(T), and the local hybridization function Hicriair) are random variables. More specifically, 
interpreting the site indices as labels enumerating the elements of random samples, Giaiaij), Sio- ('''), and Hiaia{T) 
can be understood as particular realizations of the random variables GJ^'^(t), YI'°'^{t), and H^J"^{t), respectively. Thus, 
equations (|10|1 - H16I) are in fact stochastic recursion relations, from which random samples for the random variables 
Gj^'^(r), SJ^°^(t), and H^J"^{t) can be constructed. The iterative solution of Eqs. (|1U|I -H16 |) is therefore equivalent to 
the calculation of the distribution for the random variables Gj^'^(r), SJj°'^(t), and H]^'^{t). 

Starting from an initial sample {Hi„i„{Ty\ for the local hybridization function H]^'^{t), a sample {Gio-io- (''')} for the 
local Green function G^°'^{t) is obtained by solving the quantum-mechanical many-body problem defined by S{i) [cf. 
Eqs. Hll|l - (|13|l ]. The Dyson equation (|16|l for the local Green function yields then a sample {I]i(j(T)} for the local 
self-energy Y}°'^{t), from which a new sample for the hybridization function is calculated from Eqs. H15|l . (|14|l . and 
(IIOII . Going through this loop iteratively, the random samples are successively updated until a fix point is reached. 

The iteration is numerically very involved for Bravais lattices (for instance, simple cubic, face-centered cubic, or 
body-centered cubic lattices) , in particular the construction of the lattice Green function through Eq. lfT3|) . In practice, 
the statDMFT is therefore formulated on a Bethe lattice, which is an infinite, loop-free graph with connectivity K 
(where -|- 1 is the number of next neighbours) . Besides the reduced numerical complexity, the Bethe lattice has the 
additional advantage that the statistical dependence of the random variables can be analysed in detail. 

Finally, we mention that the statDMFT reduces to the DMFT for lattices with large coordination number. For the 
particular case of a Bethe lattice, this is shown in the appendix. 



C. Specification to a Bethe lattice 



Due to the absence of closed loops, the hybridization function on the Bethe lattice depends only on the local 
cavity Green function, which can be again obtained from an appropriate effective single-site action. As a result, the 
statDMFT simplifies enormously on a Bethe lattice. 

Specifically for the K = 2 Bethe lattice with nearest neighbour coupling shown in Fig. ^ the hybridization function 
becomes 

ff»a..(i^„) = J' G%.^iLOn), (17) 

j NN i 

where now only local cavity Green functions appear. This is the crucial simplification because G^j^jg.{iujn) can 

be obtained from the effective single-site action 5^*^(j), which has the same structure as 5(j), except that the 
hybridization function is now given by 

iJ«.,(zt.„) = J2 J2 G|^';l(*'^n). (18) 
l^i NN j 
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In addition, on a Bethe lattice, because of the absence of closed loops, the exclusion of site j already ensures that site 
i is not visited, i.e., G'|^;^(iw„) = G|^]^(iw„). Accordingly, except for the 'base' site i, the hybridization function for 
all other sides is simply given by the sum over the 'forward' cavity Green function. For the base site, in contrast, all 
the neighbouring cavity Green functions contribute [see Eq. (|17|l ]. 

As a consequence, taking i as the base site, from which the cavity construction starts, gives for the Dyson equation 
for the local base site Green function 



(iuJn) 



IUJ„ 



j NN i 



(19) 



whereas the Dyson equation for the local Green function at all other sites reads 



IUJ„ 



i^i NN j 



(20) 



Note, all cavity Green functions G^*Jjg.(icj„) obey Eq. pOfl : hence they are identically distributed. Moreover, because 
the cavity self-energy 'S^^^{juj„) is obtained from 5^'^(j), containing only 'forward' sites with respect to site j, the 



cavity Green function are also statistically independent and G^j^j^^iujn) can be interpreted as a particular realization 
of the random variable G™''(ia;„). The statistical dependence between cavity Green functions at different energies 
is another issue. Naturally, it depends on the interaction processes and the particular approximation employed to 
calculate the local self-energy (see below). 

In principle the base site self-energy Eio-(iwn) and the cavity self-energy S^*J(iCi;„) are different. If, in the spirit of 

the DMFT, we neglect however this difference and identify Yiia-{i(^n) with Ej*J(i(jj„), Eq. I|19|l reduces for i and j next 
neighbour sites to 



[G. 



.PG%^{^u;^) 



(21) 



Within this approximation, the local Green function Giaiaii^n) is expressed in terms of two statistically independent 
random variables, G'fJ^^{iu!n) and G^*j^g.(ia;„). The identification of 'Eja-{i'-Un) with I]^2(*'^n) is on the Bethe lattice 
analogous to the substitution Eq. I|15|l for Bravais lattices. Even on the Bethe lattice, the statDMFT involves therefore 
two approximations: the truncation of the effective single-site action and the conflation of the interaction self-energies 
[reminiscent of the substitution Eq. (|15l) for Bravais lattices] . 



D. Localization criterion 



A natural measure of the itinerancy suggested by the effective single-site action S{i) is the total tunneling rate 
from a given site i defined by the imaginary part of the hybridization function, Ticr{u!) — js^jv^ ^ N^^J (uj), where 

Nj^J{lu) is the local cavity density of states (LCDOS), i.e., the imaginary part of the cavity Green fmiction G^j^j^{uj). 
Obviously, a finite tunneling rate rio-(w) implies an extended state at energy lu. Localized states, in contrast, lead 

to a vanishing tunneling rate. The tunneling rate vanishes if the LCDOS N^l^{ui) vanishes, which, in this sense, 
is a transport quantity, and could by itself be used as a kind of 'order parameter' for localization. The tunneling 
rate Ti^luj) as well as the LCDOS n!-'J (w) are random variables and whether they vanish or not depends on their 
distribution. The shape of the distributions has to change therefore dramatically at the localization transition. 

To understand why the shape of the distribution has to change, it is convenient to look at the distribution of a 
related quantity, the local density of states (LDOS), defined by 

N,„{lu) = --llnG^a^A^): (22) 
TT 

because it directly reflects the spatial dependence of the wave function: For an extended state at energy lu, where 
the weight of the wave functions is more or less the same on every site, the distribution of the LDOS at energy 
LU is symmetric and the most probable value coincides with the arithmetic mean value (average value). Localized 
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states with energy cj, on the other hand, have substantial weight only on a few sites. The distribution is therefore 
extremely asymmetric, with a most probable value much smaller than the arithmetic mean value, which is therefore 
not representative anymore. Since the LDOS Nja{uj) is closely related to the LCDOS N^^J{uj), the distributions 
for the LCDOS and the tunneling rate undergo the same characteristic change at the localization transition as the 
distribution for the LDOS. 

A priori it is not clear by what moment (or moments) the asymmetric distributions should be characterized. At this 
point it is useful to anticipate that distributions, in particular the distribution for the local hybridization function, 
can be used to calculate averaged transport quantities, such as the dc conductivity or the return probability (see 
below). The most probable value, i.e., the maximum of the distribution, is then expected to play a crucial role. [Note, 
averaged four-point functions (transport quantities) capture the localization effect, in contrast to averaged two-point 
functions.] It would be therefore natural to focus on the most probable values. The drawback is, however, that they 
cannot be directly obtained from the random samples; they require the explicit construction of histograms. More 
convenient quantities seem to be the so-called typical values (Dobrosavljevic and Kotliar 1997, Dobrosavljevic and 
Kotliar 1998), which are simply the geometric mean values of the random samples and yet capture the asymmetry of 
the distributions reasonably well. Therefore, we focus in this paper on the typical values. 

In particular, we consider the distribution of the LDOS and distinguish localized from extended states by a detailed 
investigation of the average LDOS, 



ivr(t^) = ^E^-H, (23) 

i 

and the typical LDOS 



exp 



(24) 



Specifically, we classify states at energy lu with N^^'^{uj) 7^ as localized if N^yP{uj) and as extended if iV*^P(ti;) 7^ 0. 
Note that the LDOS is defined for real energies. The numerical solution has to take this into account (see below). It 
should be stressed, Mirlin and Fyodorov (1994) also used the LDOS as an order parameter for the Anderson transition. 

III. NUMERICAL RESULTS 

A. StatDMFT equations for a single electron 

This section presents a detailed numerical investigation of a single electron in the AHM at zero temperature. We 
solve the statDMFT equations on a Bethe lattice where, as mentioned in the previous section, the numerical effort 
is manageable. As in other areas of statistical physics, the Bethe lattice calculation has the status of a mean field 
calculation for Bravais lattices. We expect therefore our results to be qualitatively valid for any Bravais lattice with 
d>3. 

The easiest way to obtain the zero temperature, single-particle limit of the statDMFT equations [Eqs. (|10|I -H1() () ] is 
to work with an effective single-site Hamiltonian, instead of an effective single-site action. Neglecting the spin, which 
for a single electron is irrelevant, the effective Hamiltonian reads 

HUP = ^jc]c, + nb]bj - y^ib] + b,)c]c, + E.ala, +Y,iT.Acj + H.c], (25) 

where the notation indicates that the Hamiltonian models the dynamics encoded in the cavity effective single-site 
action S{jf^ . 

The Hamiltonian representation of Sijf-"^^ is not unique. In Eq. H25I) the embedding is given by the hybridization 
with an auxiliary field described by the operator a^. The auxiliary parameters Ei, and T^j have no physical meaning, 
they just parameterize the hybridization function in terms of a spectral representation: 

H,, (zc.„) = ^ . . (26) 

For a given set of parameters {Ey, ^^I'jlj hybridization function Hjj(iiUn), the local self-energy due to electron- 
phonon coupling ^^(iwn) can be deduced from Eq. (|25|l . For finite densities, this can be done only by some approxi- 
mation, introducing further uncertainties into the approach. For a single electron, however, the interaction self-energy 
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can be obtained exactly in the form of a continued fraction (Sumi 1974, Ciuchi et al. 1997), valid over the whole 
range of parameters. Thus, for a single electron, no further approximations are required. 

To be more specific, we calculate from Ea. (|25f) the single-electron Green function for T = 0, using the continued 
fraction expansion for the electron-phonon self-energy, identify then this Green function with the zero temperature, 
single-electron limit of the cavity Green function in Eq. H2U|) . Setting iw„ -I- /i — > z = w + iry, we obtain 

G^i^)^ '—rm; ' (^7) 



[-Fii(z)] OOF 

[F,,(z-m)]-i ' 



with 



F^^{z) 



3 



(28) 



The local Green function follows from Eq. (|21|) . setting again iiOn + ij, z — uj + irj, 



Gu(z) = —J-, 7-T . (29) 

[G«(z)]-1-J2gW(z) 

According to the interpretation given in the previous section. Equations H27l) - (|29|l constitute stochastic recursion 
relations for a random sample for the local Green function Gjj°^(z). Note, Green functions with shifted energy appear 
in the continued fraction. They satisfy the same recursion relations, with energy shifted to z — qQ. 

For a fixed energy z — qQ, the local Green functions Gii{z — qCl) are identically distributed. They are however not 

independent. As a consequence, for each energy z — qQ, we have to construct a separate random sample {G^^j{z^qQ)} 
(with j = 1, Ns and q — 0, 1, M, where Ns and M are the sample size and the maximum depth of the continued 
fraction, respectively), starting from an initial sample, which we successively update via a Monte Carlo algorithm, 
similar to the schemes described by other authors (Abou-Chacra et al. 1973, Girvin and Jonson 1980), drawing 
the random variables on the rhs of Eq. H27|) from the corresponding random samples created by the iteration step 
before. Iterating this process a sufficient number of times yields, as a fix point of the stochastic recursion relations, 
a self-consistent random sample for the cavity Green function at energy z — qil. The self-consistent random samples 
are then used to directly determine average and typical values of the respective random variables, or, in the form of 
histograms, the distributions associated with them. 

For reasonable numerical accuracy, the sample size Ns , which should not be confused with the actual size of the 
Bethe lattice, but instead gives the precision with which we construct the random sample (cf. distributions), has to 
be sufficiently large. Typical sample sizes are Ng ~ 50 000. The maximum depth M of the continued fraction, which 
describes the maximum number of virtual phonons in the lattice, has to be large enough in order to capture polaron 
formation. As a rough estimate we use M w 5g^, where = Ep/il is approximately the average number of virtual 
phonons comprising the phonon cloud of the polaron. 

To distinguish localized from extended states, it is crucial to investigate the stochastic recursion relations in the 
limit ?7 ^ (Abou-Chacra et al. 1973). Numerically this seems delicate. However, initializing the iteration loop 
with a finite imaginary part of the cavity Green functions, the stochastic recursion relations can be iterated without 
problems for rj = (Alvermann 2003). The correct, self-consistent value of the imaginary part of the cavity Green 
function is then realized during the iteration. 

Tracking random samples as a function of the iteration index provides therefore information about the character 
of the states at the energy for which the sample is constructed. Random samples corresponding to extended states 
converge rapidly to its fix point, whereas random samples corresponding to localized states do not converge. Hence, 
extended and localized states can be detected by the differences they give rise to in the mathematical properties of 
the recursion relations. 

Sometimes it is advantageous to keep rj finite. But then it is necessary to monitor the fiow of the random samples, 
or of its associated average and typical values, as rj gets smaller and smaller. The part of the spectrum corresponding 
to extended states is insensitive to the ?7-scaling, after r] is small enough to see the correct LDOS nothing changes 
anymore. Localized states on the other hand are strongly affected. In fact, it can be numerically verified, that the 
typical LDOS indeed scales to zero with — > (Miller and Derrida 1993, Alvermann 2003). 

To investigate the localization properties of electron and polaron states, we specifically focused on the typical 
LDOS, tracking this quantity, for given energy and model parameters, as a function of rj and the number of iterations. 
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N. 

I 

FIG. 2: Distribution for the LDOS at a; = for 7 = 0.2,0.5, and 1.0. The inset shows the distribution for 7 = 1.5 on a 
logarithmic scale. In all cases r? = 10~^°. 

TABLE I: Average, most probable, and typical LDOS at lj = corresponding to the distributions shown in figure|5| 
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J^typ 


0.2 


0.844 


0.843 


0.841 


0.5 


0.822 


0.749 


0.773 


1.0 


0.710 


0.162 


0.445 


1.5 


0.568 


0.006 


0.106 



Introducing a rescaled transfer amplitude J = J/ Vk, we measure energies in units of the bare bandwidth Wq = 4 J = 1 
and define A — Ep/2J, a — fi/J, and 7 = 7/4J. All calculations, except in the appendix, are carried out for the 
K = 2 Bethe lattice. 



B. Electron states 

To demonstrate the feasibility of the statDMFT and the associated localization criterion, we first discuss the 
localization properties of a single electron without coupling to phonons. In the next subsection, we will then contrast 
this well understood situation with the results for the coupled electron-phonon system. Some of the results of this 
subsection are taken from Alvermann (2003). 

We begin with a qualitative discussion of the distribution for the LDOS at = 0, shown in figure El for 7 = 
0.2,0.5, 1.0, and 1.5 (inset). The shape of the distribution changes dramatically with increasing disorder: For small 
disorder (7 < 1), the distribution is approximately a Gaussian. The most probable value and the average value 
almost coincide and change little with disorder (see table P). Thus, the main effect of small disorder is to broaden 
the distribution. For larger disorder (7 > 1), on the other hand, the distribution develops a long tail, where the most 
probable value is now significantly smaller than the average value (see again table P). For even larger disorder, close 
to the localization transition (for instance, for 7 = 1.5 shown in the inset of figure[21), the distribution becomes almost 
log-normal (Montroll and Shlesinger 1983). The average LDOS is now completely meaningless. Obviously, disorder 
does not only broaden the distribution but induces, if sufficiently strong, a change in the nature of the distribution. 
From physical considerations, it is clear, that the distribution of the LDOS for a fixed energy depends strongly 
on whether the states at that energy are extended or localized. The amplitude of an extended state is uniformly 
distributed over the whole system. Thus, the LDOS does not fluctuate much from site to site, i.e., the distribution 
of the LDOS at a given energy has to be centered around an average LDOS. On the other hand, the amplitude of a 
localized state is essentially zero everywhere, except in a small range around a central site. For a given energy there 
is a certain number of degenerate localized states, each centered at a different site. The LDOS fluctuates therefore 
strongly throughout the lattice. Accordingly, the distribution of the LDOS, for a given energy, has to be very broad, 
with a long tail, due to the few sites where localized states have a large weight. Naturally, a complete characterization 
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FIG. 3: Average (upper curve) and typical LDOS (lower curve) for 7 — 1.5 and 77 = 10 The vertical dotted lines indicate 
two energies, one below and one above the lower mobility edge. 



of the LDOS has to take the changes in the distribution into account, especially near the localization transition. 

Within the statDMFT, distributions are represented by random samples. Any moment of the distribution, as well 
as the most probable value can be obtained from the histogram associated with the random sample. As indicated 
in Section II. D, in practice it is often better to avoid the construction of the histogram and to characterize, instead, 
physical quantities by their typical values, which are simply the geometric averages of the random samples. From 
Tabled we see that the typical values capture the characteristic asymmetry of the distribution in the localized regime 
reasonably well (N^yp <c N"''"^ for 7 > 1). It is therefore indeed meaningful to use the typical LDOS as a kind of 
order parameter. 

Let us now turn to a quantitative analysis of our results. In figure|31we depict, for 7 = 1.5, the typical and average 
LDOS over the whole spectral range of the pure Anderson model. The two main effects of disorder, the appearance 
of Lifshitz tails below cu — ±1.0 and the existence of mobility edges around ui ~ ±0.8, characterized by a vanishing of 
the typical LDOS (recall subsection II.D), can be clearly seen. The data have been obtained for t] = 10~^°, which is 
sufHciently small to reveal the intrinsic spectrum and to allow for a first estimate of the mobility edge. 

A precise quantitative determination of the position of the mobility edge, requires a calculation of the typical LDOS 
for ?7 — > 0. We specifically consider the lower mobility edge. With the obvious modifications all statements also hold 
for the upper mobility edge. To perform the rj ~^ limit numerically, we initialize the sample with a finite LDOS and 
then iterate the stochastic recursion relations for rj = until convergence. The iteration process yields the correct, 
self-consistent LDOS. This procedure works extremely well, as can be seen in figure^ where we plot the typical and 
average LDOS as a function of the iteration step (in other words, the sample generation) for 7 = 1.5, r] = 0, and two 
energies, one below the lower mobility edge {lo = —0.9) and one above (uj = 0). As expected, below the lower mobility 
edge, the self-consistent value of the typical LDOS continues to decrease with sample generation, whereas above the 
mobility edge, the typical LDOS remains finite. Note, in both cases, the average LDOS stays finite and is essentially 
independent of the sample generation. This is the key observation, below (above) the lower mobility edge, the typical 
LDOS vanishes (stays finite), whereas the average LDOS remains finite in both cases. 

At this point we should mention that the stable calculation of the average LDOS is subtle, because of the pure 
statistics associated with the long tail, which, on the other hand, determines the average value. To overcome this 
problem, we calculated the average LDOS after each iteration step and then performed an arithmetic average over 
the obtained values of the average LDOS. Since the typical number of iterations Nq ~ 1000, we thereby effectively 
enlarged the sample size by three orders of magnitude, which was sufficient to obtain smooth data for the average 
LDOS (Alvermann 2003). 

Thus, the tracking of the typical LDOS as a function of the sample generation can be utilized to decide whether 
states for a given energy are localized or not. In practice one would screen all energies comprising the spectrum and 
set a small threshold for the typical LDOS, for instance 10""^°, below which the typical LDOS is assumed to be de 
facto zero. According to our localization criterion, states at that energy are then classified as localized. Since this 
procedure can be performed arbitrarily close to the mobility edge, the position of the mobility edge can be determined 
very precisely. It should be also mentioned that values for the typical LDOS of the order of 10"^°" are not numerical 
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FIG. 5; The left and right panels show, respectively, the ?7-scaling of the distribution for the LDOS below {uj — —0.9) and 
above {lo = 0) the lower mobility edge for 7 — 1.5. In the insets we display the corresponding average and typical LDOS as a 
function of rj. 



artifacts. The small values are correctly captured by the floating point representation of the typical LDOS. 

An alternative is to use finite values of ij and to track the typical LDOS with decreasing ry. Below the lower mobility 
edge, the typical LDOS converges after a sufficient number of iterations to a value whose scale is set by the chosen 77 
value, whereas above the lower mobility edge, the typical LDOS reaches an intrinsic value, completely unrelated to 
rj. The average LDOS is in both cases again finite. For this procedure to work, 77 has to be of course small enough 
to resolve the intrinsic spectrum. For instance, 77 should not mimic Lifshitz tails. That the procedure based on the 
ry-dependence works excellently, can be verified in figure |S1 To clarify the origin of the differences in the Ty-scaling of 
the typical and average LDOS below and above the lower mobility edge, we first look at the distributions. Below the 
lower mobility edge at w = —0.9, shown in left panel of figure 13 the distribution changes radically with decreasing 
77 (note the logarithmic scale). The maximum of the distribution, i.e. the most probable value, shifts to very small 
values. In fact, 77 sets the scale for the most probable value. Accordingly, the typical value shown in the inset, 
decreases with 77. Note, however, the average value stays at a constant value, independent of 77. Clearly, the average 
value is determined by the long tail, i.e., the few sites, where the localized wave functions have appreciable weight. 
The average LDOS is therefore extremely less probable and indeed not a representative value. In contrast, above the 
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lower mobility edge at w = 0, shown in the right panel of figure the distribution is essentially independent of 77. 
Accordingly, the typical as well as the average LDOS do not change with rj. 

Below the mobility edge, the distributions for the LDOS, the LCDOS, and the tunneling rate r'°'^(w) become 
singular for 77 ^ 0. This behaviour is closely connected with the configuration averaged, spectrally resolved return 
probability 

Uuj) = lim(,7A)(|G,,(a; + iv)\'){e^}, (30) 

which is finite for localized states and zero for extended states. To reveal this relationship, we follow Logan and 
Wolynes (1987) and introduce a joint distribution for the real and imaginary parts of the local hybridization function, 
P(i?(w), r(w)), with R{uj) = Rei?'°=(cj) and T{uj) = ImiJ'°'=(w), in terms of which we obtain 

i^f^r^iML. ,31, 

with 

If the state at energy uj is localized, the (marginal) distribution for r(a;), i.e. P(r{uj)) — dR{uj)P{R{uj),r{uj)), is 
strongly peaked around r™P^(w) (cf. figure^. In that regime, D(w,r(w)) w S{r{uj) - r'"P"(w))iV=™(w), with 
iV=°"(w) = -(l/7r)ImGj^^(w) denoting the LCDOS (Logan and Wolynes 1987). Accordingly, Eq. (EU reduces to 

Uu;)^lim l> (33) 

V^O 1+7] ^1 "^P^(CJ) 

which, provided N''°-'"{uj) ^ 0, is finite, because 1 + i]^^r"'-P'" (uj) — > 2 for 77 ^ 0. In other words, the singular 
behaviour of the distribution for r'°'^((jj) gives rise to a pole in (16*^(0; + i'r])\'^){ei} which in turn makes fii{uj) finite. 
From diagrammatic approaches it is known that maximally crossed diagrams, responsible for destructive interference, 
produce a pole in {\Gii{u!+iri)\'^) ^^.y. Equation (|31|l provides therefore a link between the interference-based description 
of Anderson localization and the probabilistic approach adopted by the statDMFT. The ry-scaling discussed in figure[51 
can be followed all the way down to i] = 0, suggesting a precise procedure for the determination of mobility edges. In 
practice, however, we cut off the scaling at a certain threshold, for instance 10~^°. If the typical LDOS for a given 
energy cu is below the threshold, and if for that energy the average LDOS is finite, we classify the state at energy cu 
as localized. For uj = this is shown in figure |S1 Clearly, for w = the typical LDOS tends to zero for 7 ^ 2.9. The 
average LDOS on the other hand is independent of rj and remains finite. Thus, according to our localization criterion. 
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FIG. 7: Mobility edge trajectory for an electron in the pure Anderson model. The dotted lines indicate the position of the 
band edges defined by ±(M^/2 + 7/2). 



the state at energy w = is localized for 7 > 2.9, in accordance with results obtained by Girvin and Jonson (1980) 
and by Miller and Derrida (1993). 

Performing the analysis for all energies comprising the spectrum, we can systematically map out the mobility 
edge trajectory. The result is shown in figure [3 Note that the statDMFT is accurate enough to detect the two 
transitions associated with the characteristic reentrance behaviour of the mobility edge trajectory near the band 
edge: The delocalization transition at small disorder and the localization transition at large disorder. Even in the 
strongly disordered regime, where the mobility edges move to the centre of the band and all states are localized, 
the statDMFT works reliably well. We did not attempt to estimate the statistical error, but it should be of the 
order of the fluctuations visible in the plot. Our results are in excellent agreement with results obtained by Miller 
and Derrida (1993), suggesting that within the statDMFT mobility edges for noninteracting electrons can be indeed 
determined. In the next subsection, we shall demonstrate, that the statDMFT is flexible and powerful enough, to 
determine mobility edges in interacting systems as well. 



We now turn our attention to the consequences of the electron-phonon coupling. If the electron-phonon coupling is 
weak, the main effect is that the electron states above the phonon emission threshold, i.e., states whose kinetic energy 
is at least the phonon energy fi, acquire a finite lifetime due to inelastic scattering, which gives rise to lm'E,i{Lu) ^ 0. As 
a result, the critical disorder needed to locahze these states will be larger than without electron-phonon couphng. For 
the states below the phonon emission threshold, we have ImEi(a') = 0, and the critical disorder is the same as for the 
respective states in the pure Anderson model (Bronold and Fehske 2002). Increasing the electron-phonon coupling 
to the point where A > 1/VK and > 1, polaron formation starts and the LDOS fragments into an increasing 
number of polaronic subbands (Sumi 1974, Ciuchi et al. 1997). The lowest subbands are either completely or partly 
coherent, i.e. ImEi(a;) = 0, whereas the higher order subbands are strongly damped because Iml]i(a;) 7^ 0. Thus, 
in the polaron regime, we have two energy scales: the bare bandwidth and the much smaller width of the polaronic 
subbands. Obviously, the width of the distribution of the on-site energies can be small on the scale of the bare 
bandwidth but large on the scale of the subband. We distinguish therefore between two regimes: The Holstein regime, 
where disorder is small on the scale of the bare bandwidth and the Anderson regime, where it is large on the scale of 
the width of the subband. 



C. 



Polaron states 



1. Holstein regime 



For our purpose the lowest polaronic subband is of particular interest, because it is always completely coherent, 
i.e. Iml]i(aj) ~ throughout the subband, and no inelastic polaron-phonon scattering interferes with the localization 
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FIG. 8: Average and typical LDOS for the lowest polaronic subband in the anti-adiabatic, strong coupling regime for 7 = 
2.5 X W\ W = 3.45 X 10"'* is the width of the lowest subband of the pure Holstein model in units of Wo = 4J = 1. The polaron 
parameters are a = 2.25 and A = 9. The vertical dotted lines indicate two energies, one above (oj = —4.50725) and one below 
[lu — —4.50759) the upper mobility edge for the subband. 



properties of the polaron. 

First, we consider the anti-adiabatic strong coupling regime, where the phonon admixture of the states comprising 
the lowest polaronic subband is almost energy independent. Phonon- induced long-range tunneling as well as band 
flattening (Stephan 1996, Fehske and 'Wellein 1997a) are absent and the LDOS without disorder is rather symmetric. 
Numerical simulations (Stephan 1996, Fehske and Wellein 1997a) of the Holstein model indicate that the lowest 
subband in that case is just a rescaled bare band. Accordingly, we expect disorder to affect the polaronic subband 
in the same way as it affects the band of the pure Anderson model. Most notably, the mobility edges should 
simultaneously appear on both sides of the subband. Furthermore, the mobility edge trajectories should show the 
characteristic reentrance behaviour at the band edges and the critical disorder needed to localize all states of the 
subband should be determined by the states in the centre of the subband. In fact, with an appropriate scaling, the 
mobility edge trajectories for the lowest polaronic subband and the pure Anderson model should coincide. We now 
demonstrate that this is indeed the case. We start with figure |H1 which shows, for a = 2.25 and A = 9, the average 
and typical LDOS for the lowest subband; 7 — 2.5 x W with W = 3.45 x 10"^ the width of the subband in units of 
Wq = 4J. For comparison, we also show the DMFT result for the LDOS. As expected, the average LDOS is rather 
symmetric, the small asymmetry reflecting the fact that we are not yet in the extreme anti-adiabatic regime. As far as 
the effect of disorder is concerned, we see the same overall features as in the pure Anderson model (cf. figure ISJ: the 
appearance of Lifshitz tails and of mobility edges. The typical LDOS, shown in figure|Hlfor 77 = 10~* and 77 = 10~^^, 
vanishes on the both sides of the subband, clearly indicating the existence of mobility edges on both sides of the 
subband. As discussed before, a precise determination of the position of the mobility edges requires either to calculate 
the typical LDOS for = and to track the iteration flow with sample generation, or to follow the r^-scaling of the 
LDOS. For a quantitative calculation of the mobility edges, we adopted here the 77-scaling approach. In figure we 
show, for the two representative energies indicated in figure|Hl the 77-scaling of the distribution for the typical LDOS, 
the typical LDOS, and the average LDOS. Above the upper mobility edge at w = —4.50725, the distribution shows 
the characteristic properties associated with localized states: An extremely small most probable value and a long tail 
resulting in an average value much larger than the most probable value. The typical LDOS again scales to the fixed 
rj value, while the average LDOS is independent of ry. For the extended states at uj — —4.50759, on the other hand, 
we again see that the distribution is insensitive to 77 and both average and typical LDOS converge to finite values 
independent of 77. 

As for the pure Anderson model, tracking the 77-scaling of the LDOS for all energies comprising the subband, 
enables us to map out the mobility edge trajectory for the subband. From physical considerations, we know that 
the trajectory should be the same as for the pure Anderson model. Instead of calculating the whole mobility edge 
trajectory, it is therefore sufficient to verify, for representative energies, that the two mobility edge trajectories indeed 
coincide. 

As a first demonstration of this assertion, we plot in figure [TUl the typical LDOS for uj — —4.50759 (roughly the 
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FIG. 9: The left and right panels show, respectively, the 77-scaling of the distribution for the LDOS at the two energies indicated 
in figure |S] The polaron parameters are the same as in figure |H] In the left panel ui = —4.50725 (above the upper mobility 
edge), while in the right panel lu = —4.50759 (below the upper mobility edge). The insets show the corresponding average and 
typical LDOS as a function of 17. 
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FIG. 10: The typical LDOS at lo — —4.50759 (roughly the centre of the subband) scaled to its value for 7 = as a function 
of disorder measured in units of the width of the subband W; a = 2.25 and A = 9. For comparison, we also plot the scaled 
typical LDOS for the pure Anderson model as a function of scaled disorder. 



centre of the subband) scaled to its value at 7 = as a function of disorder in units of the width of the subband and 
compare it with the scaled typical LDOS at a; = for the pure Anderson model, again as a function of 7 scaled to 
the bandwidth. As anticipated, the two plots coincide. As in the pure Anderson model, the state at the centre of the 
(sub)band is localized for j/W > 2.9. To provide further evidence for our assertion that the localization properties 
of the lowest subband in the anti-adiabatic strong coupling regime are the ones of a rescaled Anderson model, we 
explicitly calculated for a few selected energies the critical disorder needed to localize the states at these energies and, 
thereby, constructed parts of the mobility edge trajectory for the lowest subband. Note the extremely high precision of 
our approach, which enables us to determine the mobility edge trajectory of a subband whose width without disorder 
is 3.45 X 10"''^ in units of Wq = 4J = 1. The results are shown in figure^2(cf. circles). Obviously, the data points 
for the subband follow exactly the trajectory of the pure Anderson model, even at very strong disorder, where all 
states of the subband are localized (see data point at j/W ~ 2.9). In the anti-adiabatic strong coupling regime, the 
localization behaviour of the states of the lowest subband is therefore identical to the behaviour of the states in a 
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FIG. 11: Parts of the mobility edge trajectories for the lowest polaronic subband of the Anderson-Holstein model in the anti- 
adiabatic strong coupling (q — 2.25 and A = 9) and the adiabatic intermediate-to-strong coupling regime {a = 0.2 and A = 1). 
Disorder is measured in units of the respective width of the subbands, which, for the plot, are scaled to Wo ~ 4J — 1. 

rescaled Anderson model. The polaron effect, in particular the band collapse, only changes the energy scale on which 
localization takes place, the internal structure of the polaron states does not affect the localization properties. Thus, 
in the anti-adiabatic strong coupling regime, disorder can be sufficiently strong to localize all states of the lowest 
subband and yet too small to interfere with the internal structure of the polaron. 

This is not the case in the adiabatic intermediate-to- strong coupling regime, where we find substantial differences 
between the localization properties of the lowest polaronic subband and the band of the pure Anderson model. 
The differences originate in the changing composite structure of the polaron states within the lowest subband. In 
particular at the top of the subband, where the band flattening due to the hybridization with the optical phonon 
branch significantly modifies the LDOS, the localization properties deviate strongly from a rescaled pure Anderson 
model. But also at the bottom of the subband significant deviations can be observed because of the phonon-induced 
long-range tunneling processes. 

In the upper panel of figure [T^ we show, for a = 0.2, A = 1, and 7 = 0.25 x W, the average and typical LDOS for 
the lowest subband. Here, W = 8.123 x 10~^ is the width of the subband for the pure Holstein model (in units of 
Wo — AJ ^ 1). In contrast to the anti-adiabatic strong coupling case, the LDOS is strongly asymmetric, a direct result 
of the band flattening which yields a steeple at the top of the LDOS. States in the steeple, belonging to the flat part of 
the dispersion, are already sluggish and therefore very susceptible to disorder. The typical LDOS vanishes therefore 
rapidly at the top of the subband. At the bottom, however, where states are rather mobile due to phonon-induced 
long-range tunneling, and hence less affected by disorder, the typical LDOS is finite. To determine the position of 
the mobility edges we again employ the ?7-scaling approach. In figure ITTI we show for a = 0.2 and A = 1 parts of the 
mobility edge trajectory for the lowest subband in the adiabatic strong coupling regime (cf. diamonds) . As expected, 
the states at the bottom of the subband are almost insensitive to small amounts of disorder, resulting in a lower 
mobility edge which is pinned at the lower band edge, which in fact does also not change. States at the top of the 
subband, in contrast, are immediately affected by disorder. Small amounts of disorder (small even on the scale of the 
subband) are sufficient to shift the upper band edge to higher energies and to localize the states at the top of the 
subband. The upper mobility edge moves very fast away from the upper band edge. As a result, the mobility edge 
trajectory for the lowest subband in the adiabatic strong coupling regime is very asymmetric (for small disorder). 

For large disorder (not shown in figure [TT|) . of the order of the width of the subband, the upper band edge of 
the disorder-broadened lowest subband moves into a spectral range with significant inelastic scattering. The reason 
is the following: Without disorder, subbands are separated by singularities in the imaginary part of the electron- 
phonon self-energy. Increasing disorder redistributes states. In particular, it moves states into the gap region, which 
separated the first from the second subband, thereby giving rise to a shrinking gap. Concomitantly, the singularity 
in the electron-phonon self-energy broadens, its state repelling character weakens, resulting in an enhanced inelastic 
scattering rate for the states at the top of the lowest subband, which strongly suppresses localization in this spectral 
range. Both effects, the merging of the subbands and the increased inelastic scattering, occur before a reentrance 
behaviour of the mobility edge trajectories can be observed. 
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FIG. 12: Average (crosses) and typical (stars) LDOS for a — 0.2, A = 1 and four different values for the disorder strength: 
7 = 0.25 xW, 7 = lxW, 7 = 2xW, and 7 = 3 x W (from top to bottom), with W = 8.123 x 10"^ the width of the lowest 
polaron subband of the pure Holstein model (in units Wo — 4 J — 1). The solid line is the DMFT result for the LCDOS. 



To support the scenario just described we depict in the three lower panels of figure ^1 for polaron parameters 
a = 0.2 and A = 1, the typical and average LDOS for 7 = 1 x VF, 2 x and 3 x W. Again, W = 8.123 x 10"^ is 
the width of the subband for the pure Holstein model (in units of Wq = 4J = 1). In all three panels, all states of the 
lowest subband are now delocalized, the typical LDOS is finite for all energies comprising the lowest subband. For the 
upper panel, where a small gap still separates the two lowest subbands, the asymmetric localization behaviour of the 
states at the bottom and the top of the lowest subband is still visible. Nevertheless, inelastic scattering is so strong as 
to make the typical LDOS finite even at the top of the subband. With increasing disorder, shown in the three lower 
panels, the gap vanishes. More importantly, however, the typical LDOS increases with increasing disorder. Thus, 
in this regime, disorder delocalizes states. Eventually, of course, with disorder on the scale of the bare bandwidth, 
all states would be localized. However, the critical disorder strength is larger than in the case of the pure Anderson 
model, because the inelastic scattering has to be overcome. 

The merging of the subbands signals that the internal structure of the polaron states and disorder start to strongly 
interfere with each other. With the vanishing subband structure, the concept of a subband mobility edge trajectory 
breaks down and a reentrance behaviour of the mobility edge trajectory cannot be established. 



2. Anderson regime 



The vanishing subband structure characterizes a transition regime, between the Holstein and the Anderson regime, 
where disorder is on the scale of the energy gap between the subbands. The localization properties in the transition 
region are rather complicated and beyond the scope of this paper. The situation becomes more transparent in what 
we call the Anderson regime, where disorder is large compared to the width of the polaronic subbands. In figure IT^ 
we show the typical and average LDOS for a = 0.2, A = 0.75, and 7 = 2. Without electron-phonon coupling, there 
would be mobility edges at w ~ ±0.9. In the presence of electron-phonon coupling, we first note that the symmetry 
between the lower and upper mobility edges is broken. The lower mobility edge is located at w « —0.8, which is 
above the mobility edge of the noninteracting system. (Here, we did not perform the 77-scaling necessary for a precise 
determination of the mobility edge.) Thus, at the low energy edge, electron-phonon coupling enhances the tendency 
towards localization. The upper mobility edge, in contrast, is shifted to higher energies (w « 1.2), i.e., at the upper 
band edge electron-phonon coupling delocalizes states and works against localization. If we mapped out a mobility 
edge trajectory, we would find a pronounced asymmetry between the lower and the upper part of the trajectory. 

The asymmetry can be explained, if we recall that our calculation is for T = 0. Therefore, states at energy lu can, 
due to electron-phonon interaction, only couple to states at energies less than lu. This leads at the high energy side 
of the LDOS to a phonon-induced coupling of localized states above the mobility edge of the noninteracting system 
to delocalized states below. As a consequence, the localized states become delocalized. On the low energy side, the 
situation is different. States below the lower mobility edge of the noninteracting system remain localized, because they 
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FIG. 13: Average (crosses) and typical LDOS (stars) in the Anderson regime for a — 0.2, A = 0.75, and 7 = 2. The solid line 
is the LCDOS obtained from DMFT (a blow-up of the low-energy part is shown in the inset). Dashed lines depict the average 
and typical LDOS for A = 0, respectively, and the vertical dotted lines indicate the mobility edges of the noninteracting system. 

can only couple to states which are already localized. Above the lower mobility edge of the noninteracting system, 
electron-phonon interaction attempts however to transform electronic band states into polaronic (sub)band states, as 
suggested by Anderson (1972). Hence, these states become heavier and more susceptible to disorder. As a consequence, 
the lower mobility edge of the interacting system shifts above the lower mobility edge of the noninteracting one. 

Further insight about the nature of the states can be gained from the LCDOS obtained within the DMFT (cL solid 
line in figure IT^ and inset). Due to electron-phonon coupling, the DMFT results for the LCDOS show at the low 
energy side pronounced plateaus with a width given by the phonon energy. (The spikes arc numerical artifacts.) The 
step-like increase of the DMFT LCDOS, together with the vanishing of the typical LDOS, is a clear signature for 
localized polaron states. The character of the states is here revealed by the DMFT LCDOS: It arises from polaronic 
defects [cf. independent boson model (Mahan 1990)], which are centered around different on-site energies, because 
of the disorder. Since they are decoupled, the DMFT LCDOS (with a higher resolution we would see the same in 
the average LDOS) does not change with energy as long as the fluctuations of the on-site energies are smaller than 
the phonon energy. If the difference in on-site energy is equal to the phonon energy, a step arises because states with 
one additional phonon contribute. The step heights reflect therefore the phonon distribution of the polaronic defect 
states. 

Whereas in the adiabatic intermediate-to-strong coupling regime disorder affects states differently at the bottom 
and the top of the polaronic subband, in the Anderson regime, electron-phonon coupling affects states in the vicinity 
of the lower and upper mobility edges of the noninteracting system differently. Only in the vicinity of the lower 
mobility edge of the noninteracting system, electron-phonon coupling and localization work in the same direction. At 
the upper mobility edge, electron-phonon coupling in fact delocalizes states. 

IV. CONCLUSIONS 

We have adopted the statDMFT to the Anderson-Holstein model to develop a self-consistent theory of localization 
in a generic disordered electron-phonon system. In particular, we investigated the localization properties of a single 
electron over a wide range of polaron parameters (adiabaticity and electron-phonon coupling strength) . 

Our approach is non-perturbative in the electron-phonon coupling and accounts for the spatial fluctuations of the 
environment by promoting all variables of the theory to random variables. The object of the theory is to calculate 
distributions (random samples). Of particular importance are the distributions for the LDOS, LCDOS, and the 
tunneling rate. Whereas information about the spatial dependence of the wave function can be extracted from 
the distribution for the LDOS, i.e. ^ P{Ni), the configuration averaged, spectrally resolved return probability is 
closely related to the distributions for the LCDOS and the tunneling rate, i.e. fa <-^- P{N^-'''), P(Ti) [cf. Ea.l31j. 

In this paper, we focused on the distribution for the LDOS. Since it is a measure of the spatial distribution of the 
wave function, it contains direct information about the localization properties: States are localized if the distribution 
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of the LDOS is very asymmetric with an extremely long tail. Although the tail, which is due to the few sites on 
which (localized) wave functions are finite, makes the average LDOS finite, the more representative typical LDOS 
vanishes. Delocalized states, on the other hand, are characterized by a symmetric distribution of the LDOS for which 
the typical and average LDOS are of the same order of magnitude. 

Using the typical LDOS to distinguish localized from extended states, we investigated in detail the localization 
properties of a single electron in the Anderson- Holstein model. We distinguished two parameter regimes: The weakly 
disordered Holstein regime, where disorder is small on the scale of the bare bandwidth and the polaronic subband 
structure is well developed, and the strongly disordered Anderson regime, with disorder large on the scale of the width 
of the subband, strong enough to interfere with the internal structure of the polaron states. 

In the weakly disordered Holstein regime, we focused on a comparison of the mobility edge trajectories for the 
lowest polaronic subband and the bare band of the Anderson model. We found significant deviations, most notably 
in the adiabatic, intermediate-to-strong coupling regime. Two main conclusions can be drawn from our studies: In 
the strong coupling, anti-adiabatic regime, the internal structure of the polaron is irrelevant and the localization 
properties of the lowest polaronic subband are essentially the ones of a rescaled Anderson model. The overall scale 
of the disorder is of course much smaller in the polaron case, suggesting that Holstein polarons are most probably 
always localized at T = 0. The internal structure of the phonon dressing contributes the most in the adiabatic 
intermediate-to-strong coupling regime. Initially, for small disorder ^ < W {W the width of the lowest subband in 
the clean system) , states at the high-energy edge of the lowest polaronic subband are extremely sensitive to disorder 
and rather easy to localize, whereas states at the bottom are almost insensitive to disorder. In contrast to the pure 
Anderson model, the mobility edge trajectory (in the disorder range where it can be defined) is not symmetric with 
respect to the band centre. Whereas the lower mobility edge is pinned to the lower band edge, the upper mobility 
edge moves quickly away from the upper band edge. A reentrance of the upper mobility edge towards the centre of 
the subband cannot be observed, because at larger disorder strength subbands vanish and the concept of a subband 
mobility edge trajectory obviously breaks down. 

The merging of polaron subbands is characteristic for a transition regime where disorder is of the order of the energy 
gap and strongly interferes with the internal structure of the polaron states. The detailed investigation of this regime 
is beyond the scope of this paper. Instead we presented results for the Anderson regime, where disorder is much larger 
than the width of the subbands. In this regime disorder and electron-phonon coupling work in the same direction 
at the bottom of the band, where electron-phonon coupling even enhances the tendency towards localization, and 
against each other at the top, where electron-phonon coupling delocalizes states above the upper mobility edge of the 
noninteracting system. In the Anderson regime, the mobility edge trajectory would acquire therefore a pronounced 
asymmetry. 

Several issues must be however clarified before real materials with polaronic excitations can be analysed along the 
lines presented in this paper. In most polaronic materials electron densities are finite. The most pressing issue is 
therefore the consideration of the polaron-polaron interaction, mediated either due to phonons or due to the Coulomb 
potential of the electron charge. In principle, the statDMFT is capable of tackling this challenging problem. The 
main obstacle here is the controlled calculation of the interaction self-energy for the ensemble of impurity models. 
The exact continued fraction expansion works only for a single electron and cannot be adopted to finite densities. For 
finite densities, the interaction sclf-encrgy has to be obtained by perturbation theory. To overcome the limitations 
of the perturbation theory, it is also conceivable to use direct numerical simulations of the ensemble of the local 
Green function, based, for instance, on Quantum-Monte-Carlo or Exact Diagonalization techniques. The localization 
properties per se can again be simply extracted from an analysis of the distribution for the LDOS. 
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APPENDIX A: DMFT LIMIT 

For a lattice with infinite coordination number, the statDMFT reduces to the DMFT. To demonstrate this point, 
it is convenient to scale the transfer amplitude J J/VK (Bronold and Fehske 2002). 

Restricting for simplicity the discussion to the Bethe lattice, Eq. shows that for if — > oo cavity and local 
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FIG. 14: Average LDOS for the Bethe lattice with connectivity K 
line is the DMFT LDOS, i.e. the LDOS for K ^ oo. 



2,4 and 8 (solid lines form top to bottom). The dashed 



Green function are identical and Eq. H19I) becomes 

Gi(jio-(iw„) = — 



(Al) 



where we have used the central limit theorem to replace the hybridization function on the rhs by the average of the 
local Green function. At this point the scaling of J is most convenient. 

The self-energy is calculated from the effective single-site action S{i) where, again due to the central limit theorem, 
the hybridization function is replaced by JG^™. The self-energy is therefore a functional of G^™(ia;„), i.e., on the rhs 
of Eq. (|Aip the only random variable is . The sample average (l/K) T,f=ii--) of Eq. ijHJ is therefore identical to 
the site average over e^. Denoting the site average by (...)£;, we get 



UUJ„ 



IUJ„ 



e,-J'G^{iUJr,)~^a[e^,G^]] ' 



(A2) 



which, together with 



S{^) 



K- 



Sint{i) 



+ e, - ix]5{t - t') + PGr(T - r') c,,(r') 



(A3) 



constitutes the DMFT equations for the average local Green function G^™(iw). For a single electron, Eqs. (|A2p 
and ljA3(l reduce, for the case without disorder, to the dynamical coherent-potential approximation (DCPA) equa- 
tions (Sumi 1974) and, for the case without electron-phonon coupling, to the coherent-potential approximation (CPA) 
equations (Elliott et al. 1974). For the pure Anderson model this is illustrated in figure mi With increasing lattice 
coordination number K we observe two effects. First, the difference between cavity and local Green function indeed 
disappears, as can be seen from the vanishing of the bumps due to the van-Hove singularities; second, the local density 
of states converges to the CPA density of states. Note also, Lifshitz tails, which are present within the statDMFT, 
also vanish, as expected, with increasing K. 
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